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Abstract 


Hawking proposed that the Sun may harbor a primordial black hole (BH) whose accretion supplies some of the 
solar luminosity. Such an object would have formed within the first 1 s after the Big Bang with the mass of a moon 
or an asteroid. These light BHs are a candidate solution to the dark matter problem, and could grow to become 
stellar-mass BHs if captured by stars. Here we compute the evolution of stars having such a BH at their center. We 
find that such objects can be surprisingly long-lived, with the lightest BHs having no influence over stellar 
evolution, while more massive ones consume the star over time to produce a range of "diga consequences. 
Models of the Sun born about a BH whose mass has since grown to approximately 107% Mo are compatible with 
current observations. In this scenario, the Sun would first dim to half its current maie over a span of 100 Myr 
as the accretion starts to generate enough energy to quench nuclear reactions. The Sun would then expand into a 
fully convective star, where it would shine luminously for potentially several gigayears with an enriched surface 
helium abundance, first as a sub-subgiant star, and later as a red straggler, before becoming a subsolar-mass BH. 
We also present results for a range of stellar masses and metallicities. The unique internal structures of stars 
harboring BHs may make it possible for asteroseismology to discover them, should they exist. We conclude with a 
list of open problems and predictions. 


Unified Astronomy Thesaurus concepts: Stellar evolution (1599); Primordial black holes (1292); Dark matter (353) 
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1. Introduction 


The dark matter problem has now become serious (Clayton 
et al. 1975). Numerous lines of evidence—such as from galaxy 
rotation curves (Zwicky 1933; Rubin & Ford 1970), the large- 
scale structure of the universe (Hernquist et al. 1996), and the 
cosmic microwave background (Planck Collaboration et al. 
2016)—indicate that most of the matter in the Universe is 
invisible. Yet despite nearly a century of research, the origin of 
this matter remains unknown, and no compelling evidence has 
emerged for a solution. Leading candidates include novel 
particles such as axion-like particles, weakly interactive 
massive particles (WIMPs), and sterile neutrinos (for reviews, 
see, e.g., Feng 2010; Tanabashi et al. 2018). 

One proposed solution is compact objects formed within the 
first second after the Big Bang: primordial black holes (PBHs; 
Carr & Hawking 1974; see Carr & Kuhnel 2021; Escrivà et al. 
2022 for reviews). These can arise from inhomogeneities in the 
initial state of the Universe and result in black holes (BHs) 
whose masses are proportional to the cube of the time of their 
formation (Carr 1975). PBHs are considered attractive because, 
unlike many other solutions, they require no modification to the 
standard model of particle physics. Carr et al. (2023) enumerate 
several hints that may point to the probable production of at 
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least some PBHs during the aftermath of the Big Bang, but it is 
currently unclear how many were produced, and whether they 
are sufficient in number to explain the dark matter. Ultimately, 
dark matter may be a combination of various solutions, 
possibly including PBHs (e.g., Eroshenko 2016; Carr et al. 
2021). 

There are several promising avenues for detecting PBHs, 
though disentangling them from more mundane astrophysical 
signals is a challenge. Gravitational-wave observations 
from LIGO/Virgo present one such opportunity (Clesse & 
García-Bellido 2017; Jedamzik 2020; Carr et al. 2023; Escrivà 
et al. 2023). Four candidate binary BH mergers may have a 
subsolar component (Phukon et al. 2021), which, if confirmed, 
could not possibly come from stars, and would be a clear 
indication of a primordial origin. There are also several with 
progenitor masses between 60 and 110M, (The LIGO 
Scientific Collaboration et al. 2021), which are within the 
pair-instability mass gap (Farmer et al. 2019; Farag et al. 2022) 
and hence cannot originate directly from the collapsing cores of 
stars. Plausible explanations include repeated mergers 
(Rodriguez et al. 2019; Fragione et al. 2022) or growth from 
smaller seeds in a nuclear star cluster (Natarajan 2021), but 
these may also be PBHs. 

Another possible route to the discovery of PBHs may be 
through gravitational microlensing. The OGLE survey has 
detected numerous compact bodies of very low mass 
(10 ^-10 ^? M.) that may be free-floating planets or PBHs 
(Niikura et al. 2019). The PBH scenario also leads to a range of 
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cosmological and dynamical consequences that may lend 
support for their production, such as cosmological background 
correlations and their effects on the reionization history 
of the Universe (Kashlinsky et al. 2005, 2018; Cappelluti 
et al. 2013, 2022; Kashlinsky 2016; Boldrini et al. 2020; 
Hasinger 2020; Carr et al. 2023). 

The solar system also provides some tools for detecting 
PBHs. The Earth, Moon, and Sun can be used as transient PBH 
detectors, as typically assumed velocities would have PBHs 
passing through these objects faster than escape velocity, which 
would lead to unique dynamics (Li et al. 2023), oscillations 
(Kesden & Hanasoge 2011; Luo et al. 2012), and craters 
(Yalinewich & Caplan 2021). Scholtz & Unwin (2020) have 
argued that the hypothesized Planet 9 (Batygin & Brown 2016) 
could be a PBH, which could be confirmed by annihilation 
signals. 

The milky Way is expected to contain ~100 million BHs 
from normal stellar evolution pathways, with the average BH 
being at a distance of 21 pc (7-106 au) (Sweeney et al. 2022). If 
PBHs at the classical Hawking evaporation limit (~10~'* Mo) 
constitute the dark matter, this number increases to ~10°° BHs 
at an average distance of ~l au. This would be far more 
numerous and far more densely spaced than stars, raising the 
possibility of their capture by stars (Ilie et al. 2021; Ilie & 
Paulin 2022) or star-forming clouds. 

If PBHs comprise all or part of the galactic dark matter halo, 
then they may have orbits with random inclinations and 
velocities from a Maxwell-Boltzmann distribution. PBHs in 
the slow tail of this distribution, if aligned with the orbits of 
stars or star-forming regions, have some chance of being 
captured. Montero-Camacho et al. (2019) have found that the 
probabilities of capture by a star are low, as a PBH on average 
falls in faster than the stellar escape velocity and hence transits 
the star with only negligible dissipation from dynamical 
friction and accretion. The capture of PBHs during star 
formation is much more likely due to the time-dependent 
gravitational potential of the adiabatically collapsing cloud 
(Capela et al. 2013, 2014; Eroshenko 2023). PBH capture in a 
star-forming region is even more likely in dwarf galaxies due to 
their lower mean velocities as well as in the early Universe, 
which has formed the basis for additional observational 
constraints on PBH dark matter (Oncins et al. 2022; Esser & 
Tinyakov 2023). 

A low-velocity PBH captured by a star would slowly 
consume it (Clayton et al. 1975; Picchio 1981; Oncins et al. 
2022). At early times, the PBH may accrete at a Bondi-like rate 
(dMgyu/dt ox Mj, with Mpy the mass of the BH), where 
accretion is limited by the speed of sound in the center of the 
star. It may later transition to an Eddington-like accretion 
(dMgy/dt x Mgy) if the radiation pressure can stall the freefall 
accretion, but this is uncertain. 

Following Hawking's suggestion that the Sun may harbor a 
PBH (Hawking 1971), and inspired by the then-unresolved 
solar neutrino problem (Stothers & Ezer 1973), Clayton et al. 
(1975) created the first solar BH models from partial 
evolutionary sequences. Thorne et al. (1981) and Flammang 
(1982, 1984) considered optically thick, Bondi-type accretion 
and the role of gas pressure, and Marković (1995a, 1995b) 
analyzed the effects of convection, rotation, and magnetic 
fields. On the other end of the mass spectrum, BHs inside of 
supermassive stars (M > 1000 Mo), quasi-stars, have been 
evoked to explain the origin of supermassive BHs in the early 
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Universe (Begelman et al. 2008; Begelman 2010; Volonteri & 
Begelman 2010; Ball et al. 2011, 2012). Several authors have 
studied the growth of PBHs inside neutron stars (Kouvaris & 
Tinyakov 2014; Baumgarte & Shapiro 2021; Richards et al. 
2021; Schnauck et al. 2021). 

In this work, we carry these models forward and consider the 
evolution of stars having 0.8-100 Mo that have formed about a 
BH spanning all possible masses up to the stellar mass, 
including the first full numerical evolutionary simulations of 
solar-like stars with a central BH. We utilize a 1D stellar 
evolution code that makes use of modern opacities and other 
microphysics in the calculations, as well as multiple treatments 
of the accretion. We describe the physics and implementation 
in Section 2. We present detailed numerical results on solar 
models in Section 3.1, and analytically explore a grid of stellar 
models in Section 3.2. We discuss the limitations of our model 
in Section 4, and list several open problems. Finally, we 
conclude the paper in Section 5. 


2. Methods 


In order to simulate the evolution of a star under the 
hypothesis that it has a central BH, we treat the BH as a point 
mass at the core and solve the 1D stellar structure equations 
with modified inner boundary conditions (Clayton et al. 1975; 
Ball 2012; Farmer et al. 2023) through an extension to the 
MESA (Modules for Experiments in Stellar Astrophysics, 
123.05.1; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn 
et al. 2022) stellar evolution code. 

In the standard equations of stellar structure, the mass M, 
radius R, and luminosity L are all zero at the center point. We 
set the mass at the inner boundary to Mgy plus the mass of a 
cavity around the BH from which it accretes. This models the 
Bondi sphere, where the infall velocity is greater than the speed 
of sound in the medium (Bondi 1952). We take the Bondi 
radius as the radial inner boundary condition, defined as 


GM 
Rg = 2 ZE, (1) 


Cs 


where c, is the adiabatic speed of sound at Rg. We take the 
mass in this region as (87/3) pRg, assuming a density profile of 
por ?? (Ball 2012). 

The luminosity generated from accretion is given by 


€ dMpu E 
l—e dt 


L= 


: (2) 


where c is the speed of light and e is the radiative efficiency. As 
a first case, we study a fixed radiative efficiency at the 
canonical €— 0.08 roughly corresponding to reaching the 
innermost stable orbit around a Schwarzschild BH, but this 
efficiency is highly uncertain (see Section 4.4). For compar- 
ison, the efficiency of nuclear fusion is an order of magnitude 
smaller at ~0.007, and the efficiency of a rotating Kerr BH can 
be much larger. 

We take the luminosity as the lesser of the Eddington 
luminosity Lg and a Bondi-like convective luminosity Lg 
(Ball 2012), defined as 


Lg = 4r E GMpu, (3) 
K 
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Lg = 1675 a (GMpp)?, (4) 


stl 


where & is the opacity and I’, is the first adiabatic exponent. The 
efficiency of convection is denoted 7 and we take it to be 0.1 
(Ball 2012). We also consider a case without Eddington-limited 
accretion, in which the luminosity is only given by Equation (4). 
We adopt the widely used OPAL (Opacity Project at Livermore; 
Iglesias & Rogers 1996; Rogers & Nayfonov 2002) equation of 
state and opacities. To give typical values, standard models of 
the present Sun have p—150gcm ^, c,—5x 10’cms |, 
&— 15cm? g !, and T1 5 /3. Note that the metal-rich plasma 
of the solar core is much more opaque than the commonly 
adopted  electron-scattering opacity for ionized hydrogen 
(& — 0.4 cm? g’. Values for zero-age main-sequence (ZAMS) 
stars across our chosen mass range are given in Appendix A. 

The mass of the system decreases over time by the 
equivalent of the energy converted into radiation. We limit 
the time step to Mgy/(dMgy/dr) in order to stably resolve the 
accretion. We include a small amount of exponentially 
decaying convective overshoot in order to smooth convective 
boundaries. We do not consider element diffusion or mass loss 
by stellar winds. We also do not model the late stages of the 
consumption of the star. Our implementation and models are 
publicly available"? (Bellinger 2023). 


3. Results 
3.1. Solar Models 


The Sun is formidable and destroying it is generally considered 
to be a difficult problem. Stellar evolution models without a 
central BH suggest that the Sun will evolve into a black dwarf and 
live on essentially forever, possibly only destroyed by proton 
decay in the far distant future (Laughlin 2007; Caplan 2020). As a 
fiducial first case, we characterize the evolution of the Sun with a 
central PBH until its premature demise. 

We calibrated solar models with a central with masses above 
the classical Hawking evaporation limit (Mpu > 10 79 Mo). These 
models have the observed present-day properties of the Sun, but 
get some fraction of their luminosity from accretion onto a BH. 
For each PBH mass, we adjusted the initial helium abundance and 
mixing-length parameter until the models reached 1 Ro and 1 Lo 
with [Fe/H] — 0 within a tolerance of 107 at the solar age of 
4.572 Gyr. Models with an initial Mgy ; 10^ !! Mo converged; 
models with greater mass are either over-luminous or consume the 
Sun before its present age. The converged models correspond to 
PBHs that formed in the first femtoseconds of cosmic time: 
approximately 10 2^—10 ^ s after the Big Bang, following the 
inflationary epoch and preceding the quark epoch. 

We find that a PBH spanning a wide range of masses (a 
current Mg < 10? M.) could exist inside of the present Sun 
without significant modification to its present structure. 
Because the masses of these BHs are so small and the 
luminosity from their accretion is much less than the nuclear 
luminosity, they change the speed of sound, density, and 
nuclear reaction rates in the Sun by much less than a percent. 
These models are therefore compatible with present-day 
helioseismic and solar neutrino observations. This mass 
window coincides with the PBH mass range that could still 
constitute the entirety of the dark matter (Carr & Kuhnel 2021). 
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Figure 1 shows the evolution and fate of the Sun if it had 
formed about a PBH with an initial mass of 10 !! Mo. Over the 
lifetime of the Sun from the pre-main-sequence (PMS) until the 
present day (4.572 Gyr), the PBH increases its mass tenfold as 
it accretes the solar plasma. The BH drives convection in the 
core of the Sun and mixes its innermost regions, but otherwise 
causes little change to its outward appearance. The fate of the 
Sun, on the other hand, changes dramatically. By the age of 
7 Gyr, the BH consumes 0.146 of the solar mass. The solar core 
now cools, causing nuclear reactions to cease. Soon thereafter, 
the luminosity from accretion starts to vastly exceed what was 
once the luminosity from nuclear fusion (see Figure 2). 

Next, the large accretion luminosity causes the star to 
become fully convective and puff into a red giant. The partially 
fused core gets fully mixed back into the envelope, leading to a 
star with an extremely high surface helium abundance. The 
expansion of the Sun halts at a maximum of ~0.03 au, and thus 
Earth is saved from being engulfed by the red giant Sun. 
Earth's oceans nevertheless still boil off as its blackbody 
temperature reaches upward of 530 K (—250?C). 

As acoustic stellar oscillation frequencies are sensitive to the 
speed of sound of the core and hence the core temperature, the 
initial transition to lower core temperatures would likely bear 
an observable seismic signature (e.g., Unno et al. 1989; Aerts 
et al. 2010; Basu & Chaplin 2017), thus providing an early 
warning for the changes that are to come. The discovery of 
solar g-modes would also help to falsify this model, as these as 
yet undetected oscillations would grant insights into the deepest 
layers of the Sun. 

In this model, the Sun has about 8 Gyr after the present day 
remaining before it will become a BH, including a >2 Gyr time 
span in which the BH mass is greater than 196 of the solar 
mass. In the Kippenhahn diagram (Figure 1), two distinct 
changes of slope are evident in the BH growth: at 5 Gyr 
(Mag > jg-9 M) when the accretion transitions from the 
Bondi rate to Eddington; and again at 7 Gyr (Mpu ~ 10 ^ Mo), 
when the accretion takes over to supply the full luminosity of 
the star. 

Figure 2 additionally shows the evolution of the Sun if it 
formed about a PBH with a greater initial mass of 107" Mo. 
Here the situation is quite different, as it would have quenched 
nuclear reactions 2 Gyr ago and already expanded into a 
giant star. 

Figure 3 shows the evolution of these models in the 
Hertzsprung-Russell (H-R) diagram. The unique and long- 
lived post-main-sequence (MS) phase (i.e., after fusion has 
ceased) may motivate observational searches. In normal MS 
evolution, the luminosity grows and the star migrates up the 
H-R diagram. These models follow a rather different path. The 
accretion energy causes the star to evolve essentially backward 
along the PMS track of standard single-star evolution, though 
persisting in those stages for much longer. As nuclear reactions 
begin to shut off, these stars initially migrate toward lower 
luminosities and cool to temperatures as low as 4300 K. After 
the complete cessation of nuclear fusion, the star then grows in 
luminosity and radius, while remaining up to hundreds of 
kelvins cooler than the normal red giant branch (RGB). The 
temperature in this phase depends on the composition, which in 
turn depends on how much helium was produced during the 
evolution of the MS before entering the accretion-driven phase. 

Stars in this sparsely populated region of the H-R diagram 
are known as sub-subgiants and red stragglers. Hundreds of 
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Figure 1. Kippenhahn diagrams showing the evolution of the interior of the Sun with and without a central BH. The left panels show the mass distribution, with 
regions of energy generation and transport indicated. The right panels show the radial distribution, with the radius of the photosphere (black line) and the solar radius 
(horizontal dashed line) indicated. The top panels correspond to a normal solar evolution model evolved through the MS until core hydrogen exhaustion and up 
through hydrogen shell burning as a red giant. The bottom panels show a model that is consistent with the present Sun with a BH growing at its center. Nuclear fusion 
(red) provides the bulk of the solar luminosity until the BH is of sufficient mass to quench the reactions. The BH drives convection (hatches), which mixes the 
innermost parts of the core, and eventually the entire star. Note the differences in the y-axis scale between the panels. 


examples are now known from detailed studies of color- 
magnitude diagrams (Geller et al. 2017; Leiner et al. 
2017, 2022) and many are known to be RS CVn binaries. 
Leiner et al. (2017) created evolutionary models showing three 
plausible explanations for these stars: binary mass transfer, 
envelope stripping, and magnetic activity. Nevertheless, they 
form an interesting sample in which to potentially search for 
stars harboring BHs. 

The evolution of a star with a central BH in the red straggler 
phase may proceed slowly over multiple gigayears, with the 
timescale set by the accretion onto the BH. In contrast, normal 
red giants climb the RGB proceeds over a span of less than 
1 Gyr, with evolution driven by the contracting core and the 
presence of a hydrogen-burning shell (Miller Bertolami 2022). 
The models eventually exceed ~10 Lo and ~5 Ro. 

Due to the fully convective nature of the post-MS, 
candidates may be identified from spectra of RGB stars with 
large He abundances. The Sun was about 70% hydrogen at its 


birth, and will convert ~10% of its mass to helium over a 
period of 10 Gyr. If the Sun became fully mixed after the MS, 
its surface helium abundance would rise to about 35%. This 
material otherwise normally stays trapped in the core; low-mass 
stars also generally have a lower surface helium abundance 
than at their birth due to element diffusion and gravitational 
settling (e.g., Christensen-Dalsgaard et al. 1993). 

Additionally, such stars may also potentially be identified by 
an absence of mixed modes in their asteroseismic oscillation 
spectrum. This is because g-modes, and hence the mixed 
modes that are typical of red giants (e.g. Hekker & 
Christensen-Dalsgaard 2017), require stable stratification that 
is ordinarily provided by the dense radiative helium core. An 
accounting of such modes, as for instance has been done with 
depressed dipole modes (Fuller et al. 2015), may already be 
useful in constraining capture rates. 

Finally, Figure 3 also shows that for sufficiently low BH 
masses, the star evolves through these phases indistinguishably 
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Figure 2. MESA simulations of the past and future luminosity of the Sun assuming no central BH (left), a central BH of initial mass 10~'' Mo (center), and one of 
10? Mo (right). The solar age and luminosity are indicated as well as the mass of the BH. Unlike the first two models, the model on the right is obviously not 
compatible with the present Sun. In both of the right panels, it takes ~1 Gyr for the BH to transition from supplying 1% to supplying 100% of the solar flux. When the 
BH approaches the point of supplying 20% of the solar luminosity, nuclear reactions begin to wane and eventually cease. The Sun then dims to half its current 
luminosity over the span of 100 Myr. Over the course of the next 300 Myr, the BH grows to a mass of 5 x 10? Mo, at which point it produces the full solar 
luminosity. Soon thereafter, the Sun expands into a red giant, and does so at a much greater pace than in the normal case: instead of taking 6 Gyr for the Sun to double 
in luminosity, it is achieved with the aid of the BH in ~300 Myr. Finally, the accretion luminosity drives the Sun to slowly expand over the course of several gigayears 
until it produces a maximum of more than 10 Lo before the simulation terminates. The BH in the bottom of the middle panel having a mass of 10 ? Mọ is depicted 
with its actual Schwarzschild radius (~3 mm); more massive BHs are too large to show (e.g., ~30 m at 10 ? Mo) on this plot. 


from a star without a BH inside it. For solar mass stars with 
Mpgy = 107 i Mo, the star deviates from normal solar evolution 
only at the base of the RGB; for even smaller masses, the star 
survives the ascent up the RGB without significant BH growth. 
Future studies will be required to understand the post-main- 
sequence evolution of these objects, such as if they still ignite 
the helium flash and eventually evolve toward the white dwarf 
stage (see also the discussion in Section 4.5). 

The models stop running when the BH reaches a mass of 
approximately 0.05 Mo, at which point the Bondi radius is 
more than half the stellar radius and growing rapidly. The 
system has converted about half a percent of its total mass into 
radiation through the accretion process. Ball et al. (2011, 2012) 
found an upper limit on the ratio of the inner BH mass to the 
total mass of 0.119 for supermassive quasi-stars, which they 
related to the Schónberg-Chandrasekhar limit. We may only 
speculate as to the evolution of the system beyond this point; 
we will return to this question in Section 4.5. 

Finally, we consider a second case in which the luminosity is 
not limited by the Eddington limit, but rather by the maximum 
flux that can be carried away by convection. In this case, we 
find that the overall evolutionary path is the same as in the first 
case, only vastly accelerated. Rather than the long ~8 Gyr 
BH-driven evolutionary phase, the star puffs into a red 
straggler and subsequently consumes the star over a time span 
of less than 100 Myr. After the BH reaches a critical mass, the 
Sun reacts to the large injection of central energy by expanding 
into a giant over a thermal timescale. The evolution then 
continues over thermal timescales, which are in turn con- 
tinuously shortened due to the rapidly increasing surface 
luminosity. Rather than the evolution stopping at ~15 Lo as in 
the Eddington-limited case, here the star reaches 10* Lo before 
the simulation ends. Another difference from the Eddington- 
limited case is the fate of Earth: while in the previous case, it 
was somewhat uncertain, in this case, Earth is almost certain to 


be engulfed, as the maximum radius of the expanding Sun 
grows to exceed an astronomical unit. As in the case of normal 
solar evolution, these details may change depending on the 
efficiency of wind-driven mass loss, which is neglected in these 
simulations. 


3.2. Stellar Models 


We have seen in the previous sections that solar models can 
be surprisingly long-lived for a wide range of PBH seed 
masses. Low-mass BHs have no effect on the evolution of the 
Sun, whereas more massive BHs bring about a unique, 
accretion-driven evolutionary phase. If the accretion luminosity 
is limited by the Eddington luminosity, then this phase may last 
several gigayears; otherwise, the star is consumed in tens of 
megayears. 

In this section, we address lifetimes more generally across a 
range of stellar masses. Because the pre-MS phase is relatively 
short and accretion onto the BH during this phase is slow, we 
focus here on the MS. Similar to the analyses by Roncadelli 
et al. (2009) and Oncins et al. (2022), we have computed a grid 
of ZAMS models, and solved Equations (3) and (4) analytically 
to determine how long it takes for a BH at the center of these 
models to fully accrete the star. The details of the model grid 
and calculations can be found in Appendix A. Here we focus 
on a discussion of the main findings. 

We must first emphasize that unlike in the previous sections, 
which presented detailed numerical simulations for 1 M. 
evolutionary tracks, these analytic estimates ignore the 
evolution of the internal structure of the star, as well as any 
interplay between the accretion luminosity and the stellar 
plasma. These approximations therefore become particularly 
questionable when the luminosity from the accretion exceeds 
the luminosity from fusion. Therefore, numerical simulations 
will be necessary to confirm these statements for additional 
masses other than the solar mass. That being said, we have 
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Figure 3. H-R diagram showing the MS and subsequent evolution of the Sun with and without a central BH. The solid black line representing the normal case 
becomes dotted when the Sun has exhausted its central hydrogen abundance. The tracks with a BH become dotted when the accretion luminosity exceeds the 
luminosity from fusion. In all cases, the star expands into a red giant. Stars with lower BH masses than those shown proceed through these evolutionary phases as 
normal. Ages and the ZAMS are indicated; notice that a star with a BH dims rapidly before climbing redward of the RGB slowly over 7 Gyr. The normal Sun climbs 


the RGB much quicker over a span of about 50 Myr. 


found that the numerical results in Section 3.1 for the 1 Mo 
case are in excellent agreement with the results found below, 
even for large BH masses. For these estimates, we assume the 
growth of the BH is limited by the Eddington luminosity at a 
fixed radiative efficiency as in Section 3.1. For the second case, 
the results would be the same at early times, but the lifetimes 
would be cut short around the time the Eddington limit is 
reached. 

Figure 4 shows the lifetimes of stars harboring a central BH. 
The mass at which the Eddington limit is reached spans from 
about 10~'° Mo for subsolar mass stars to 10 7 Mo for 20 Mo 
stars. According to the first accretion scheme, the MS lifetimes 
of massive stars (M > 10 Mọ) would not be cut short at all 
unless they are born about a BH with a mass greater than 
~15% of their mass, even when the BH accretion is supplying 
the full stellar luminosity. A star's longevity relative to the 
mass of the BH increases with increasing stellar mass. These 
estimates give that a 20 Mo star born about a 9 Mo BH would 
live as long as the normal MS lifetime of a 20 Mo star. Low- 
mass stars are especially long-lived; a 0.7 Mo star formed about 
a BH having 1% of its mass may live ~1Gyr. The 
characteristic shape of the figure is dictated by the transition 
from Bondi to Eddington accretion, explaining for example 
why the 10 and 1 Gyr contour lines have opposite concavity. 


This figure furthermore shows us which kinds of stars are 
sensitive to which masses of PBHs. MS stars with very low 
mass (0.2 < Mą/Mo < 3) are potentially sensitive probes of 
very low-mass PBHs do P < Mpgu/Mo < 107°). Solar mass 
stars for example would be good probes of ^10 !? Mo PBHs 
because such objects reduce the expected lifetime of the star by 
half if captured during formation, thus giving ample time to 
witness a Sun in the midst of being devoured. On the other 
hand, high-mass stars have very low ratios of density to the 
speed of sound (i.e., low erodibility), and therefore BHs in their 
center are unable to accrete efficiently in the Bondi regime. 
Massive MS stars (M, > 10 Mo) are thus only sensitive to 
more massive BHs (Mpu > 1 Mo) over the short lifetime of 
their MS. However, as mentioned, these analytic estimates 
ignore post-MS evolution as well as the detailed interaction of 
the BH with the stellar plasma. The opacity in the core drops 
significantly as it becomes degenerate, leading to potentially 
much more rapid accretion in the Eddington regime. These 
detailed considerations will be the subject of future studies. 


4. Open Problems 


Future work is needed to develop the theory of stellar 
evolution with central BHs if this method is to show promise of 
detecting or constraining PBHs in the asteroid-mass window. 
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Figure 4. Analytic accretion estimates for the lifetimes of MS stars born harboring a PBH. Contours of constant age are given by black lines; for low-mass BHs, the lifetimes 
of the stars are unchanged (brown shading), but as the BH mass approaches the stellar mass, the ages go to zero (blue shading). The colors indicate the percentage of the 
normal MS lifetime that the star will live before being swallowed by the BH; brown on the top indicates that the star will survive at least as long as its ordinary MS lifetime, 
and white indicates that the BH cuts the star’s lifetime short by half. The thick white-dashed line shows the Bondi—Eddington transition, i.e., when the BH is massive enough 
for its accretion luminosity to reach the Eddington limit. The white lines in the middle show a comparison between the luminosity from accretion and the luminosity from 
fusion, showing for example that a 1 Mo star could attain 1 Lo from slowly falling into a BH of 10 * Mo. The white-dashed lines on the right show when the BH is 1%, 
10%, and 50% of the stellar mass, showing for example that a 1 Mo star born about a BH of 0.1 Mo could live for 2300 Myr. 


Here we list a few open problems that we consider to be well- 
posed and well-motivated. How many BHs were made by the 
Big Bang? How often do they become bound to star-forming 
clouds? Once captured by a cloud, how long does it take for 
them to reach the center of a star? How much energy is 
generated by their accretion of the stellar material, and can this 
luminosity be used to power the star? What are the final stages 
of stellar consumption? And what would be the consequences 
if this scenario turns out to be common? We will now 
contemplate these open problems in some detail. 


4.1. What is the Initial Mass Function of PBHs? 


Throughout this work, we have considered the idea of PBHs 
as being the solution to the dark matter problem. Even if they 
are not the full solution, they may nevertheless exist in smaller 
numbers, as natural scenarios appear to produce them in the 
early stages after the Big Bang. 

Many works that study PBHs use an idealized mass function of 
PBHs, in which it is assumed that all of those compact bodies are 
of the same mass; i.e., a so-called monochromatic mass function. 
In realistic scenarios, however, non-monochromaticity is unavoid- 
able. Every known PBH production mechanism yields a spectrum 
of above-threshold overdensities, which then collapse to form 
BHs over a range of masses (see Section II of Escrivà et al. 2023 
for an extensive discussion of various formation scenarios). For 
the most-studied case of PBH formation from inflationary 
overdensities, the phenomenon of critical collapse (Choptuik 1993; 


see also Evans & Coleman 1994; Koike et al. 1995; Niemeyer & 
Jedamzik 1998; Kühnel et al. 2016 for later work) inevitably 
broadens any original mass spectrum, even if the initial spectrum 
was monochromatic. 

After their formation, which usually happens at very early 
times, such as around the end of the quark epoch around 10 $ s 
after the Big Bang, the PBHs also begin to accrete. Since these 
are distributed throughout the Universe, being locally inhomo- 
geneous, the increase in mass will further broaden the initial 
mass and spin spectrum. The currently most minimal and 
natural scenario uses critical Higgs inflation (Garcia-Bellido & 
Ruiz Morales 2017; Ezquiaga et al. 2018) in combination with 
the thermal history of the Universe. This leads to an extended 
PBH mass function with multiple bumps at scales of the order 
1075, 1, 10, and 105M.., corresponding to the electroweak and 
quantum chromodynamic scales, the pion plateau as well as to 
e'e annihilation, respectively, while still predicting a modest 
DM fraction (710^?) in the asteroid-mass window (Carr et al. 
2021). This mechanism has been proposed to explain the 
observations mentioned in the Introduction, and could 
potentially account for the entirety of the dark matter (see 
Figure 38 of Carr et al. 2023). 


4.2. What is the Capture Rate of PBHs in Stars and Star- 
forming Clouds? 


The capture rate of PBHs by stars is likely to be quite low 
(Montero-Camacho et al. 2019) From simple kinematic 


THE ASTROPHYSICAL JOURNAL, 959:113 (14pp), 2023 December 20 


arguments, a PBH in the Galactic halo falling onto a star from 
infinity should be accelerated above the escape velocity of that 
star, quickly pass through the star, and escape. The damping 
forces from accretion and dynamical friction of a small PBH 
are so weak that there is an effectively negligible chance that 
the PBH will be captured. Only in a neutron star are these 
damping forces great enough to produce an appreciable capture 
rate, which forms the basis of constraints on PBHs from 
neutron star survivability in the Galaxy (e.g., Baumgarte & 
Shapiro 2021). 

In contrast, capture during star formation is made possible by 
the time-dependent gravitational potential of the collapsing 
cloud. During collapse, a PBH present in the cloud is not 
accelerated by the gas around it, and so it can become bound to 
the cloud (Capela et al. 2013, 2014). This capture mechanism 
depends on the PBH velocity distribution, with only those 
PBHs in the slow end of the tail being captured. This generally 
favors capture in lower-mass dwarf galaxies with smaller 
velocity distributions over those in the Galactic disk (Esser & 
Tinyakov 2023). One might also expect a higher rate of 
captures for lower PBH masses, given the greater possible 
number density. It is also possible that PBHs are clustered (see, 
e.g., Meszaros 1975; Carr 1977; Carr & Silk 1983; Freese et al. 
1983; Chisholm 2006, 2011; García-Bellido & Clesse 2018; 
Trashorras et al. 2021 or Section II of Carr et al. 2023 for an 
extensive discussion) which may allow a star to capture 
multiple PBHs. If transport toward the core is fast, it is not 
inconceivable that clustered PBHs may even merge in the 
stellar center. One might also expect a higher rate of captures 
for lower PBH masses, given the greater possible number 
density, making stars a powerful tool for probing the most 
elusive PBHs having the lowest masses. 


4.3. Can a PBH Reach the Stellar Core within Stellar 
Lifetimes? 


We have assumed that a PBH is already at the center by the 
time the star forms. However, unless the PBH was already 
located at the center during the initial collapse of a gas cloud, the 
PBH may have to undergo an in-spiral phase from the outer 
envelope to the stellar core. There are two main sources of drag 
that drive the orbit to decay: dynamical friction and hydro- 
dynamic drag. If the orbiter accretes mass, this would also play a 
role as a drag force of roughly Mv, where v is the speed of the 
orbiter. However, if M is the Bondi rate, the overall gravitational 
drag would be described by the expression of dynamical friction 
(Cantó et al. 2011; Lee & Stahler 2011, 2014). Because 
hydrodynamic drag is proportional to the cross section of the 
orbiter, this contribution is subdominant in small objects like 
PBHs. Even if the effective interacting surface area is 
comparable to nRj, the hydrodynamic drag may be at most 
comparable to the dynamical friction. Thus, for an order of 
magnitude estimate of the orbital decay timescale, we will 
assume dynamical friction to be the main source of decay to the 
orbit of the PBH. 

To take a concrete example, for a PBH with a mass of 
Mpg = 10 1? Mo circularly plowing through an envelope layer 
with a local density of p = 1074 g cm ? at a distance r from the 
center of a 1 Mo MS star, the orbit decay timescale at that 
location due to dynamical friction (Ostriker 1999) may be 
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estimated as 


InAY!( M n 
tss = 100 My s) rear | 


(Ms V(r Vp) Y os 
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where Men is the enclosed mass at r. Here A is defined as the ratio 
between the effective linear size of the surrounding medium and 
that of the perturbing object. We take the former distance scale to 
be rœ 1 Ro and the latter to be Rg. With this choice, In À very 
weakly depends on Mgu, r, and Men. For this particular case, this 
estimate suggests that the decay time would be short compared to 
the life of the star once the PBH is reasonably deep inside the 
envelope. However, it is important to remark that this estimate 
only gives an instantaneous decay time based on the local 
properties of the star. The total decay time would depend on a 
number of different factors, such as radiation from the accretion 
onto the perturber (Park & Bogdanović 2017), equation of state 
(Khajenabi & Dib 2012), magnetic field (Sánchez-Salcedo 2012), 
orbital motion (Kim & Kim 2007), and density gradient (Sánchez- 
Salcedo & Brandenburg 2001). 

That being said, the dependence of taecay on Mpy and p also 
suggests that low-mass PBHs, or those PBHs that could not 
penetrate sufficiently deep inside the envelope, could spend an 
exceedingly long time orbiting outside the stellar core— 
possibly longer than the lifetime of the star. This possibly long 
orbit decay time could have interesting implications. If the orbit 
decay time is not sufficiently shorter than the lifetime of the 
star, the evolution of the star would be essentially the same as 
that without a PBH at the core. In that case, the evolution of the 
star would in turn determine the fate of the PBH. The PBH's 
orbit decay would be significantly affected by the changes in 
the local density and speed of sound as the star evolves. If the 
star undergoes an explosive event (such as a supernova) while 
the PBH orbits outside the core, the PBH could be ejected 
along with unbound ejecta. Although the wandering PBHs may 
not affect the evolution of the star, its nonlinear orbit motion 
would generate a trailing pressure wave and overtake its own 
wake (Kim & Kim 2007), which could affect the oscillations in 
the envelope. It is, however, difficult to speculate beyond this 
point; detailed hydrodynamical simulations would be a route 
toward making further progress on this topic. 


4.4. What is the Growth Rate and Radiative Efficiency of a 
Microscopic BH? 


In this work, we have considered two accretion schemes, 
both using a fixed radiative efficiency. A Bondi-like accretion 
rate is assumed at early times. In the first scheme, the 
]uminosity is eventually limited by the Eddington luminosity; 
in the second, the luminosity is limited to the maximum that 
can be carried by convection, and thus the Bondi-like accretion 
rate is maintained. We will now discuss some of these 
considerations. 

Bondi-type accretion is possible when the specific angular 
momentum at the Bondi radius is smaller than that at the 
innermost stable circular orbit (ISCO). Assuming rigid rotation, 
the typical angular velocity of the Sun's interior is roughly 
3 x 10$ s^!. The BH mass below which the specific angular 
momentum at the Bondi radius (Rgj Q, where Q is the orbital 


THE ASTROPHYSICAL JOURNAL, 959:113 (14pp), 2023 December 20 


frequency) is smaller than that at ISCO (/12 GMgy/c) can be 
estimated as 


C 4 Q -1 
Mpu & 1 Mo s . 6 
as of 5 x 107 cm =) (; x 10-6 =) (6) 


where c, is the speed of sound at the Bondi radius. This justifies 
Bondi-type accretion at early times. 

As the flow transits toward the Eddington limit, the higher mass 
fluxes could slow down the outward diffusion of the radiated 
photons, thus trapping them and advecting along with the 
accretion flow. If the timescale for photon diffusion is larger than 
the freefall timescale at the Bondi radius, then the photon trapping 
may suppress the radiated luminosity to below ~Zg, and the gas 
can continue to accrete beyond Lg/c? (Begelman 1978). As an 
order-of-magnitude estimate, the time for photons generated in the 
vicinity of the PBH to reach the Bondi radius can be estimated via 
the photon diffusion time taie ^" Rp7/c, where Tœ pgRgk is the 
optical depth and pg is the density at the Bondi radius. This can be 
compared with the freefall time tayn at r œ Rg via 


fait || 3 Mag K 
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This order unity ratio indicates that the photon trapping may be 
important in determining the radiative efficiency, especially for 
a massive PBH embedded in a dense core. Therefore, limiting 
the mass accretion rate to the Eddington rate may be unrealistic, 
and the BH may accrete the rest of the star much faster as in the 
second case that we considered. 

There are yet further difficulties with assuming the 
Eddington limit. Flammang (1982, 1984) found that when 
gas pressure P, dominates over radiation pressure P,, as is the 
case in the solar interior, the Eddington luminosity may be 


lowered to 
Lp = dfi : zi is (8) 


1/*g8 


which at solar conditions is approximately 0.196 of the Eddington 
luminosity. If this limit holds, this implies that a BH of the same 
mass has much less influence on the surrounding plasma, and so 
the lifetime of any BH-driven phase of evolution would be 
considerably briefer. However, this reduced luminosity ignores 
the influence of rotation. Even the relatively slow rotation period 
of the Sun implies a nearly maximally spinning central BH once it 
has accreted a significant amount of mass. The angular 
momentum of the infalling material spins up the BH, slowing 
the growth of the BH and raising the luminosity above Lp 
(Marković 1995a, 1995b). It is conceivable that a disk or torus 
may form, enabling even greater luminosity. However, the 
turbulent convection resulting from the accretion luminosity 
could prevent or disrupt the formation of such structures 
(Marković 1995a). A significantly larger luminosity than the 
Flammang luminosity could also be directed along the axis of 
rotation, and a shear-induced instability could induce turbulent 
heating and further slow the accretion rate (Marković 19952). In 
summary, the theory of spherically symmetric accretion suggests 
very high accretion rates with very low accretion luminosities are 
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possible, but the complex interplay of rotation, convection, 
magnetic fields, and the potential for disk or torus formation 
complicates the picture significantly. 

If the BH does radiate at the Flammang luminosity and 
grows at the Bondi rate, then nearly all the infalling material 
solely contributes to the mass growth of the BH and almost 
none is converted into energy, and thus, there is nearly no 
pushback against the infalling stellar material. The BH still 
grows very slowly at early times, but upon hitting a critical 
mass, grows exponentially to consume the star extremely 
quickly. Since the luminosity would be low, the stellar 
properties would be unaffected until the BH is very massive, 
at which point the star has little time left before being 
consumed (Marković 1995a). Hence, in this model, the Sun 
would not enter a long-lived accretion-powered red giant phase 
if the BH reached this critical mass on the MS. 

On the other hand, limiting the luminosity at the inner 
boundary to either the Eddington or Flammang luminosity is 
questionable. Sanyal et al. (2015) studied hydrodynamical 
simulations of massive stars, and found that luminosities can 
locally exceed the Eddington limit by a factor of a few without 
driving outflows, suggesting that even much larger luminosities 
inside the star are possible. In models of supermassive 
quasi-stars, Begelman et al. (2008) take the limiting luminosity 
to be the Eddington limit of the entire system, rather than of the 
BH alone. These considerations motivated our second accretion 
scheme, in which the luminosity is limited only by the maximum 
flux that can be carried away by convection. 

Lastly is the question of the radiative efficiency itself, which 
is highly uncertain, and depends on the properties of both the 
inflowing plasma and the BH itself. For instance, the radiative 
efficiency is expected to undergo variations on the orders of 
unity in relation to the spin of the BH (Bardeen et al. 1972). As 
the mass of the BH grows by orders of magnitude, if the 
accretion coherently adds angular momentum to the BH, then 
the BH spin would change rapidly, possibly approaching a 
limiting value, beyond which any further change may not have 
much effect. 

It is unclear at present how an asymmetrical flow or 
magnetic fields may affect the accretion. The differentially 
rotating envelope convection zone of the Sun spins up a 
magnetic dynamo (Spruit 2002; Charbonneau 2020), which 
complicates the situation even further. Marković (1995a) 
considered how such magnetic fields may affect the formation 
of an accretion torus. Ultimately, more detailed simulations that 
can probe the conditions inside the Bondi sphere will be 
required in order to understand these complex phenomena, 
including the growth rates, radiative efficiencies, and super- 
Eddington accretion. 


4.5. How Does a PBH Affect the Post-MS, and What Sort of 
Transient is Associated with the Final Stages of Stellar 
Cannibalism by a Central PBH? 


As a star becomes a white dwarf, the core density increases 
along with the core electron conductivity. This drives the 
central opacities very low, which in principle would enable 
rapid accretion of the degenerate matter onto the BH 
(Roncadelli et al. 2009). However, this ignores the feedback 
from the stellar material, which may regulate this process and 
prevent rapid accretion. We plan to address the later stages of 
stellar evolution with detailed numerical simulations in 
future work. 
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Our simulations end when the Bondi radius is halfway to the 
photospheric radius, corresponding to the BH having eaten 
about 5% of the stellar mass. The fate of the remaining material 
is therefore outside the domain probed by our simulations, and 
we may only speculate beyond this point. It is unclear whether 
the BH simply accretes the rest of the star (Ball et al. 2011). 
Instead, disk formation and outflowing jets could be possible, 
which would potentially carry material and angular momentum 
away from the BH. The object may also ultimately evolve 
toward an X-ray transient and have properties similar to the 
observed population of accreting X-ray binaries. If that is the 
case, then hydrodynamic simulations of the transient phase 
may be able to identify characteristics that could distinguish 
them from those in binaries. Another potential observable is the 
disappearance of a giant star, for which upcoming transient 
searches such as LSST may be sensitive. 

The final BH spin is highly uncertain and beyond the scope 
of the simulations presented here. Nevertheless, there are many 
interesting possibilities. As an order of magnitude estimate, if a 
PBH consumes a | Mo star rotating with a period of 1 month, 
the final Kerr spin parameter a = J/(c Mgy) (with J being the 
spin angular momentum) is of order 1 km, comparable to the 
Schwarzschild radius. This suggests that PBHs that consume a 
majority of their host star could be near maximally spinning. 
On the other hand, simulations of core-collapse supernova have 
shown that lower rotation rates may arise if there is significant 
mass loss (Fuller & Ma 2019) such as if some of the angular 
momentum of the infalling gas is lost through winds 
(Belczynski et al. 2020) or jets (Gottlieb et al. 2023). Another 
consideration is if the PBH initially rotates rapidly. While most 
past work argues for small PBH spins at formation (e.g., De 
Luca et al. 2019), some calculations predict that PBHs near the 
classical evaporation limit may be spun up by emitting 
Hawking radiation, suggesting that the smallest PBHs may 
have large spins initially (Calzà & Rosa 2022). 


4.6. What Would Be the Consequences if This Scenario Turned 
Out to Be Common? 


In the relatively unlikely scenario that PBHs are both 
common and frequently captured by stars or star-forming 
regions, then they may contribute to several open problems. 

Li et al. (2023) recently found that the number of low-mass 
stars increases with metallicity. Low-mass BHs consume low- 
mass stars much more efficiently at low metallicity (see 
Figure 6) because their cores are much less opaque, thus 
increasing the Eddington limit, and therefore could contribute 
to this effect. 

Some globular clusters show stars with helium abundances 
of up to 40%, which are difficult to explain with current models 
(e.g., Lee et al. 2005) Studies of white dwarfs in globular 
clusters also indicate that they have evolved from helium-rich 
progenitor stars (Althaus et al. 2017). If early stars harboring a 
central BH drive a stellar wind after mixing the helium core 
into the envelope, then they would enrich the interstellar 
medium with helium for future generations of star formation. 

There appears to be a number of missing white dwarfs 
relative to expected numbers in some star clusters like the 
Hyades (Tinsley 1974; van den Heuvel 1975; Weidemann et al. 
1992; Portegies Zwart et al. 2001; Tremblay et al. 2012). If the 
accretion of degenerate plasma is much more efficient due to its 
very low opacity, then a large population of very low-mass 
BHs could reveal itself in this way. The Hyades also appears to 
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host several BHs (Torniamenti et al. 2023), the presence of 
which may require either very small supernova kicks or PBHs. 

The LIGO/Virgo collaboration has recently discovered 
numerous BH mergers, with some residing in the pair- 
instability mass gap (Abbott et al. 2020a, 2020b). PBHs could 
in the first place naturally originate with masses within the pair- 
instability gap. Alternatively, these BHs could originate from 
PBHs of significantly lower mass that subsequently populated 
stars and accreted them from within. 

If PBHs constitute the dark matter, then likely PBH 
production scenarios may be capable of explaining most of 
the LIGO/Virgo observations without a stellar origin. How- 
ever, a subpopulation of mergers with both components having 
^10 M. appear to still require a different explanation than 
PBHs (see Figure 29 of Carr et al. 2023). While these of course 
most likely originate from stellar channels, the mechanism 
discussed in this work could have additional bearing on this 
matter: by converting early stars into BHs, thereby increasing 
their mass, a part of the pronounced peak around a solar mass 
could shift toward larger masses, hence potentially yielding 
additional mergers in the unexplained range. 

It is possible that there is a component of the Galactic dark 
matter that co-orbits with the disk, a so-called dark disk. Some 
particle theories of dark matter with weak self-interaction allow 
for the necessary dissipation to form a disk and would be one 
component of a larger dark sector (Fan et al. 2013). Stellar 
capture of PBHs can convert PBH dark matter from the halo 
into stellar-mass BHs in the disk of the Milk Way, providing a 
natural formation mechanism for a Milky Way dark disk. 
Present constraints from Gaia on thin dark disks (scale heights 
h S; 50 pc) suggest as much as 1% of the Galactic dark matter 
could be in the disk (Buch et al. 2019; Widmark et al. 2021). A 
constraint on the PBH capture rate in the Milky Way could 
potentially be obtained from astrometric measurements. 

Recent JWST observations show over-luminous objects at early 
times (Liu & Bromm 2022; Ilie et al. 2023). Normally, a 1 Mo star 
must wait 12 Gyr to reach 10 Lo; a BH with an initial mass of 
10^? Mo could drive it to that luminosity in less than half the 
time. The dark star scenario has been evoked to explain these 
observations (Spolyar et al. 2008; Rindler-Daller et al. 2015), 
which are hypothetical stars containing WIMPs or self-interacting 
dark matter such as neutralino dark matter. On that note, we also 
wish to draw attention to numerous other works that have 
implemented candidate dark matter solutions into stellar evolution 
codes. Vincent et al. (2015) studied asymmetric dark matter in 
solar models, Martins et al. (2017) extended these to other solar- 
type stars, and Rato et al. (2021) performed a similar study in 
subgiant stars. Battich et al. (2016) and Córsico et al. (2016) 
studied axion cooling in pulsating white dwarfs, and Severino & 
Lopes (2023) carried out a similar study on the red supergiant 
Betelgeuse. Lopes et al. (2019) looked at the effects of WIMPS in 
red clump stars, and argued that a differential analysis of nearby 
giants and giants in the Galactic Center could potentially be used to 
distinguish stars harboring dark matter. Ayala et al. (2020) studied 
dark photons in RGB stars. More research will be required to 
determine which observational predictions are unique to each 
theory. 


5. Conclusions 


Here we have numerically computed the first full solar 
evolution tracks with a central BH, and also presented analytic 
results for a wide range of stellar masses. Based on our adopted 
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accretion schemes, we find that stars harboring BHs can have 
surprising longevity, which broadly shows that stars are 
compatible with an enormous population of small BHs 
inhabiting the galaxy. Every MS star living now could have 
formed about a PBH of less than 10? M, without changing 
the observed population. The present Sun could currently 
harbor a PBH as massive as the planet Mercury which would 
elude present detection capabilities. There is, however, no 
positive evidence that a BH is indeed currently in the Sun. If 
the Sun did form about a PBH, it could not have been more 
massive than about 107" Mo, depending on the adopted 
accretion scheme and radiative efficiency, as it would otherwise 
change the observed properties of the present Sun. 

As stars harboring BHs are long-lived— possibly even when 
the BH constitutes a significant fraction of the mass and 
supplies a large luminosity—this opens a number of possible 
avenues for searching for the presence of such objects. One 
promising approach is through asteroseismology. Stars like the 
Sun are gardens of sound, ringing with acoustic oscillations 
with characteristic periods of ~5 minutes. The frequencies of 
these oscillations depend intimately on both the global and 
interior structure of the star, and can be used for example to 
determine the age of the star (e.g., Christensen-Dalsgaard 1984; 
Bellinger et al. 2016; Hon et al. 2020). 

When the BH mass is small, the star is essentially 
indistinguishable from a normal star. However, the accretion 
onto the BH causes the core of the star to become convective 
(see Figure 1) mixing the material and causing the mean 
molecular weight gradient in that region to become flat. The 
cores of low-mass stars are ordinarily radiative, and therefore 
are predicted by stellar evolution theory to have steep mean 
molecular weight gradients. Asteroseismology is sensitive to 
the shape of this gradient (Cunha & Metcalfe 2007; Beck et al. 
2012; Deheuvels et al. 2015, 2016; Gehan et al. 2018) and has 
been used to characterize the masses of convective cores 
(Angelou et al. 2020). Asteroseismology has provided 
inferences into the near-core structure of some solar-type stars 
(Bellinger et al. 2017), including one with a convective core 
(Bellinger et al. 2019). A low-mass solar-type star with a 
convective core could be a signature of a central BH, as this is 
not produced by ordinary stellar evolution pathways, and could 
potentially be detectable via asteroseismology. 

A further opportunity for the discovery of a star harboring a 
PBH becomes possible if the accretion luminosity can serve as 
the star's dominant source of energy. Here the star first 
becomes a sub-subgiant and later a very cool red giant, known 
as a red straggler. Ordinary red giants have a compact radiative 
core encased in a convective envelope and oscillate in so-called 
mixed modes, which is a coupling between pressure mode 
oscillations in the envelope and gravity mode oscillations in the 
core (e.g., Hekker & Christensen-Dalsgaard 2017). However, 
these latter kinds of oscillations require stably stratified regions 
somewhere in the star; if the star is fully convective, then no 
gravity modes and hence no mixed modes are possible. 
Therefore, a giant star pulsating in pure pressure modes may 
also be a signature of a star harboring a BH at its center. 

On the other hand, if the accretion is radiatively inefficient, then 
the star would experience almost no outward change to its 
appearance until relatively soon before it is destroyed 
(Marković 19952). The star would appear as normal until at 
most thousands of years before its destruction, at which point its 
luminosity would begin to increase precipitously. These objects 
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could then potentially be discovered by their rapid disappearance. 
It is also conceivable that the stellar envelope could be ejected in 
the end stages, leading to a subsolar mass BH in a nebula. 

In a future paper, we aim to perform a detailed asteroseismic 
characterization of stars being powered by PBHs. If they 
present a unique signature, then these objects could potentially 
be discovered through the data archives of the CoRoT 
(Auvergne et al. 2009), Kepler (Borucki et al. 2010), and 
TESS (Ricker et al. 2015) missions. Currently, there is high- 
quality data for approximately 100 solar-type stars and 
thousands of giants; presently available data may already be 
of sufficient quality to find such an object. By the end of the 
decade, orders of magnitude more asteroseismic targets will 
become available thanks to the upcoming PLATO (Rauer et al. 
2014; Miglio et al. 2017) and Roman (Gould et al. 2015) 
missions, including asteroseismology of ~1 million red giant 
stars near the Galactic Center. This presents an opportunity to 
either discover such objects, or to place bounds on their number 
and capture rate. The implications for stars in more advanced 
evolutionary stages, numerical results for stars of different 
masses and metallicities, and investigations into stellar 
populations will also be explored in future works. 

There is a question of what to call these hypothetical low- 
mass quasi-stars, should they ever be found. They might be 
called “Hawking stars.” 
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Appendix A 
ZAMS Models 


To facilitate accretion timescale estimates, we show the 
central speed of sound and central opacity of our models at 
ZAMS in Figure 5. Bondi accretion is inversely related to the 
speed of sound and Eddington accretion is inversely related to 
opacity. The central speed of sound only changes by a factor 
~2 with increasing stellar mass around a typical value of 
10° cms_!. The opacity decreases with increasing stellar mass, 
eventually asymptoting for M» 10M. at a minimum of 
0.3 cm? B For comparison, past authors have assumed fully 
ionized H using & —or/m, = 0.4 cm’ g~', with Thompson 
cross section c and proton mass m,. Central densities vary by 
about two orders of magnitude across stellar masses. 

We find these central values are reasonably well fit by 


c, = 1.6 x 107 (3.5 + log; Mx), (A1) 
k —2034-12x Mz, (A2) 
—19 
logo p = F2, (A3) 


1 + exp{—2.7(log,)(Mx) — 1)} 


where M, is in Mo units and the rest are in cgs. Metallicity has 
at most a 15% effect on the speed of sound, whereas a high 
metallicity can more than double the opacity in the core of a 


THE ASTROPHYSICAL JOURNAL, 959:113 (14pp), 2023 December 20 


14 

g 

E 1.2 

9 D 
2 = 
= E 
10 9, 
E 2 
o S 
Bos e 
El 

o 

[7] 


= 
o 


1 10 100 “4 


M/Mo 


M/Mo 


Bellinger et al. 


100 


10 


density [g/cm?] 


10 100 1 10 
M/Mo 


100 


Figure 5. Central speed of sound (left), opacity (center), and density (right) for ordinary ZAMS models as a function of mass and shown for 1%, 10%, 100%, and 


150% of the solar metallicity. 


l 
I 
[L| 1 
Bondi 
Eddington 
Transition 
Lo 


accretion luminosity [Lo] 
Leda — Leonai [erg/s] 


10 5 410 7 10° 10° 
black hole mass [Mo] 


10° 


Eddington 


10° 
black hole mass [Mc] 


-11 
10 0.1 1 10 


100 
stellar mass [Mo] 


Figure 6. The transition between Bondi and Eddington accretion onto the central BH for ZAMS. Left panel: the accretion luminosity of a 1 Mo star at ZAMS (white-dashed 
line) as a function of the mass of the PBH that it formed about. Middle panel: the difference in Bondi and Eddington luminosities at ZAMS. Right panel: the transition mass 
between Bondi and Eddington accretion at ZAMS. For small BHs, the accretion proceeds at the Bondi rate (Equation (4)); larger BHs accrete at the Eddington limit 
(Equation (3)). More massive stars have significantly higher Bondi—Eddington transition masses. Metallicity is very important in determining this transition in low-mass stars. 


low-mass star. Note that we have neglected the Z dependency 
in developing these formulas. 

We use these formulas to estimate the transition point between 
Eddington and Bondi accretion for a ZAMS star by equating 
Equations (3) and (4). The results are shown in Figure 6. This 
relation at solar metallicity is well approximated by 


Mggr = 4 x 1070M}. (A4) 


All stars that form about a PBH whose mass at capture is at 
least 1—2 orders of magnitude lower than the Mger will survive 
at least the MS; and stars above 3 Mo survive even if the BH 
mass is significantly greater than the Mggr. We aim to assess its 
post-MS fate in future work, including the interesting feedback 
and self-regulation processes arising between the interplay of 
the accretion onto the BH and the response of the stellar 
plasma. Note that the densities of stellar cores increase 
substantially throughout the evolution of the MS; for instance, 
the solar core had a ZAMS density of ^80 gem >, whereas 
now it has a density of ~150 g cm ^. We ignore these effects in 
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this analytical estimate, but they are incorporated in the 
numerical one presented in the main text. 


Appendix B 
Analytic Timescale 


We are now prepared to estimate the time it takes for a PBH 
to fully accrete the star, and compare this to the MS age. To 
reasonable accuracy, MS ages follow 

tus = 10 Mz? Gyr. (B1) 


For the Bondi case, we will estimate the time it takes tgonai for 
the BH mass to grow until it reaches the Eddington limit via 


Meer 1 
J Mà 


Magno BH 


l l E ÍBondi 
Meger — Mano SBondi 


where Mgr is the BH seed mass, Mger is the mass of the 
Bondi-Eddington transition, and the Bondi accretion mass 
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sensitivity is given by 


T) 
Sgondi = — — A, (B3) 
(1 — e)n 161G^p 
For the Eddington case, we have 
Mx 
fo amu = HH, (B4) 
M, Meu Stud 


where M,, = max(Mpn,o. Mpgr) and the Eddington accretion 
mass sensitivity is 


€ CK 
Sead = ; B5 
me 1 € 4nG eu 
which implies that 
M. 
trad = Sgaa log (55) (B6) 


Finally, we used these analytic estimates to construct Figure 4. 
The characteristic shape of the curve in this figure (i.e., the 
point at which the curve begins to transition from red to blue) is 
given by equating the MS lifetime tms with the total accretion 
lifetime. 
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